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^ , ' Abstract 
O ' 

We consider the single-particle velocity distribution ol a one-dimensional fluid ol inelastic particles. Both the freely evolving 
(cooling) system and the non-equilibrium stationary state obtained in the presence of random forcing are investigated, and 
special emphasis is paid to the small inelasticity limit. The results are obtained from analytical arguments applied to the Boltz- 
mann equation along with three complementary numerical techniques (Molecular Dynamics, Direct Monte Carlo Simulation 
Methods and iterative solutions of integro-differential kinetic equations). For the freely cooling fluid, we investigate in detail 
' the scaling properties of the bimodal velocity distribution emerging close to elasticity and calculate the scaling function asso- 
ciated with the distribution function. In the heated steady state, we flnd that, depending on the inelasticity, the distribution 
function may display two different stretched exponential tails at large velocities. The inelasticity dependence of the crossover 
velocity is determined and it is found that the extremely high velocity tail may not be observable at "experimentally relevant" 
inelasticities. 
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I. INTRODUCTION 
A. Motivation 
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C ' In the widely studied context of nonequilibrium stationary states, granular gases stand out as an interesting model 
system, accessible to and subject of many experimental and analytic investigations. Their theoretical description 
and understanding is one of the presently important issues of the development of the out-of-equilibrium statistical 
mechanics. 

The main difference between molecular gases and granular gases stems from the fact that at each collision between 
e.g. steel or glass beads (in experiments), or idealized smooth hard spheres (in analytical and numerical investigations), 
a fraction of the relative kinetic energy is lost [|l|. This inelasticity is responsible for many interesting phenomena, 
such as the appearance of spatial heterogeneities, or non-Gaussian velocity distributions... Theoretically, two opposite 
situations have been extensively studied in the context of smooth inelastic hard spheres we shall consider here. Namely, 
the free cooling case where no forcing mechanism compensates the energy loss due to dissipative collisions (see e.g. 
the review and references therein), and the uniformly heated system where an external random force acts as a 
■ heating process on the grains, allowing a non-equilibrium stationary state to be reached [^|-|5|. 

In this work, we will study the two above situations (i.e. with or without energy input), and concentrate on a 
one-dimensional granular fluid. For the homogeneously heated gas (section |l|), the focus will be on the high energy 
tail of the velocity distribution P{v). Whereas the velocities up to the thermal scale obey a Maxwell-Boltzmann-like 
. distribution, we will show combining kinetic theory arguments and numerical simulations (both Molecular Dynamics 
^ ' and Monte Carlo) that in the limit of vanishing inelasticity, P{v) displays a exp(— w'^) large v behaviour. At finite 
Q inelasticity, this tail is asymptotically dominated by the law exp(— w"^/^) already predicted in M. These predictions will 
O also be confirmed by the results of a high precision iterative solution of the non-linear Boltzmann equation. On the 



other hand, without energy injection (section III), we will similarly concentrate on the limit of small inelasticity (that 
appears quite singular, unlike in the heated case), and shed some light on the importance of spatial heterogeneities 
and velocity correlations: detailed scaling properties of the solutions of the homogeneous Boltzmann equation will be 
obtained analytically and checked numerically. Further confrontation against Molecular Dynamics simulations will 
show that the velocity distributions of the Boltzmann homogeneous cooling state share some common features with 
those obtained by integrating the exact equations of motion. 



B. The model 



We shall consider a one-dimensional gas of equal mass particles of length a and density n, evolving on a line of 
length L = N/n with periodic boundary conditions. These particles undergo binary collisions with conservation of 
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momentum but loss of a fraction (a^) of the kinetic energy in the center of mass frame: consequently, if vi and V2 
(resp. v'l and v'2) are the velocities of the two particles involved before (resp. after) the collision, then 

v'l + V2 = Vl + V2 (1) 

v[-v'2 = -a{vi~V2), (2) 

where < a < 1 is the restitution coefficient. We also introduce the inelasticity parameter e = 1 — a (e = for elastic 
collisions) . 

We will focus on the behaviour of the velocity distribution P{v, t) in two cases: 

• without energy injected, the above collision rules define a system where energy dissipation through collisions 
is not balanced and the typical velocities of particles progressively decrease. This free cooling regime has 
been widely studied [^-^ and in particular in dimension 1 by molecular dynamics studies |^J^,^,^. Slight 



modifications of the collision rule allow to bypass the inelastic collapse 20 1 and to observe an asymptotic scaling 
regime for P{v, t) jlH . 

• a steady state can be reached if the loss of energy through collisions is balanced by an injection that can be 
achieved through a random force rj{t) acting on each particle 

^=i^ + 77(t), mv(t'))^2D5{t-t') (3) 

where D is the amplitude of the injected power and F the systematic force due to inelastic collisions. Velocities 
consequently execute a random walk in-between the collisions and in the coUisionless case P{v, t) obeys a 
diffusion equation with a "diffusion" coefficient D. This model was first introduced and discussed by Williams 
and MacKintosh ^ in dimension 1, and studied in higher dimensions H]; variants have also been proposed 

MM- 

We define the granular temperature as the average kinetic energy of the system: 

T{t) = / dvv^P{v,t) . (4) 



The function T{t) decreases for the freely cooling system, but it eventually fiuctuates around a steady-state value in 
the heated case. 



C. Investigation methods 



Our study relies on three complementary approaches: 

• Molecular Dynamics (MD) simulations p3| integrate the exact equation of motion in a finite box: we consider 
TV hard rods of length ct, on a line of linear size L, with periodic boundary conditions, random initial velocities, 
and we use an event-driven algorithm to study their dynamics. 

• the Boltzmann equation describes the evolution of the one-particle distribution function P(v,t), upon the 
molecular chaos factorization hypothesis [p4| . This equation is therefore a mean-field approximation of the 
problem. It can be solved numerically by the Direct Simulation Monte Carlo (DSMC) method |^^, or in certain 
cases with an even better precision by an iterative method similar to that used in |26[| . 

• in the elastic limit a ^ 1, an analytical scaling approach can be used to study the Boltzmann equation. It is 
important to note that, in the particular case of one-dimensional hard spheres, the elastic case e = is quite 
peculiar. Indeed, for e = 0, the collisions only exchange the velocities of the particles: this system is therefore 
unable to thermalize and is equivalent to one with transparent particles where P{v, t) is frozen in time. 
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II. STEADY STATE OF THE HEATED FLUID 



The exact solution of the problem of a d = 1 inelastic gas appears inaccessible, which prompted us to carry out 
molecular dynamics (MD) simulations to calculate P(v, t). In order to see whether a mean-field type approach can give 
a reasonable description of the inelastic gas, we reconsider the kinetic theory of the process. We solve the appropriate 
Boltzmann equation by simulations in the general a case, and derive exact results in the £ — > limit. 

Since the gas systematically reaches a stationary state with a temperature T that depends on the inelasticity (see 
below), it is useful to introduce the rescaled velocity 
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and the corresponding distribution function 



f{c) = ^/T{r)P{v,t). (6) 



In order to compare the velocity distributions at various inelasticities, numerical data for the rescaled velocity distri- 
butions will always be displayed in the figures. 



A. General a case 



1. Molecular Dynamics simulations 



The molecular dynamics simulations are carried out with = 5000 hard rods, using an event-driven algorithm, 
and submitting the rods to random kicks at a frequency that remains much higher than the collision frequency, in 
order to simulate the noise of equation (|^). Note that, since the Langevin equation considered in j2l] is different from 
ours 0, comparing the two approaches is not possible without further investigations on either side. 

Starting from an initial distribution of velocities, a steady-state is reached after a transient, and velocity distributions 
can be measured. Figure Q displays such distributions for inelasticity ranging from a — 0.6 up to 0.999. Strong non- 
Gaussian behaviour is observed, with over- or under-populated high energy tails depending on a. Moreover, the 
inset shows that the system remains quite homogeneous at low inelasticity, while strong density fluctuations develop 
at larger inelasticity. The value of the density n of the system seems to have no influence on the shape of the 
velocity distributions, density fluctuations seem to increase roughly proportional to l/{na) (at constant a). Detailed 
investigations of the spatial correlations are left for future studies. 
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FIG. 1. Velocities distributions for MD simulations with 5000 particles, density na — 0.5, for restitution coefficients 0.6, 0.9, 
0.95, 0.99, 0.999 (from top to bottom). The symbols show the Gaussian distribution. In the inset are shown the space density 
fluctuations for 0.99 (continuous line) and 0.6 (dotted line). 



*It contains a friction term in velocity space whose amplitude is linked to the force amplitude, whereas detailed balance is 
not satisfied by our model, where the forcing is independent of the state of the system 
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2. Kinetic theory 



Assuming that the density of the particles is low and neglecting both velocity and spatial correlations of colliding 
partners, the following Boltzmann equation for the spatially averaged velocity distribution function P(v,t) can be 
written as 

oo oc 

dtP-DdlP = -n J dv'\v - v'\P{v)P{v') + j^^^ J dv'\v-v'\P{v')P ( ^"^ ■ (7) 

— oo — oo 

The right-hand side above contains the collision terms corresponding to the "dissipative" rules , while the Fokker- 
Planck term Dd^P on the left-hand side takes into account the energy injected by the random forces (^). 

The system described by Eq. (|^) is expected to relax to a steady state since the power input is independent of the 
velocities while the loss of energy is roughly proportional to the energy itself. This expectation can be made more 
quantitative by deducing the equation for the temperature (T = (f^)) of the system 

^^2D^^{l-a'){\v-vf) (8) 

where (...) denotes averaging with P(w, t), and |w — u'| represents the relative velocity of two randomly chosen particles. 
There is a stationary solution to this equation that has a simple physical meaning. Namely the rate of input of energy 
D) is equal to the rate of loss of kinetic energy in the center-of-mass frame ~ {1 — a'^){v — v'Y (with the extra 
factor n|w — ii'l coming from the collision rate). One may also estimate the characteristic time of reaching the steady 
state. Indeed, the quantities T^^"^ = (w^)'^/^ and {\v\^) are expected to have the same leading large-time {t — > oo) 
dependence, and thus, up to an unknown numerical constant C, Eq. (^) can be written as 

dT 

— ^2D -C n(l- a^)T^''^ (9) 
dt ^ ' ^ ' 

The typical relaxation time is then Troiax ~ [D/{n(l — a^))]^/'^. In the small inelasticity limit (e 0), this relaxation 
time diverges as Troiax ~ £~^/^. We have indeed observed such a behaviour of the approach to the steady state both 
in MD and DSMC simulations. 

Many pieces of information on the stationary distribution function have been obtained in B; in particular, the 
deviations from a Gaussian $(c) = e~'^ /^/•\/27r have been investigated by the Sonine expansion 

/(c) = $(c) |^l + f;ap5p(c2)j , (10) 

where the Bp's are polynomials orthogonal with the Gaussian weight $. The coefficients Op are then obtained from 
the moments of /. From the definition of temperature, ai vanishes and the first correction 02, which is related to 
the kurtosis of the velocity distribution, has been computed in any dimension neglecting non-linear contributions of 
O (02); in dimension 1, it has the expression 

_ 4(z;4) ^ 16(1 -2a^) 
^ ~ 3(t.2)2 129 + 30a2 ' ' 

We note a peculiarity of dimension 1: 02 does not vanish as a — > 1, unlike in space dimensions c? > 1 in which 
limQ,^! 02 = 0. This is a hint that the quasi-elastic limit is more singular in c? = 1 than in higher dimensions. 

Besides it was shown in [|| how to determine the high energy tail of the velocity distribution. We briefly recall the 
argument. The e-dependent gain term in the collision integral appearing in the Boltzmann equation (|^) is a priori 
neglected. In the large velocity limit the resulting equation for the steady-state distribution Ps{v) reads 

^ -n\v\Ps{v) (12) 

which yields a high energy tail of the form exp(— |-\/n7^|wp^^). Then one verifies that the gain term is indeed 
negligible (as would be the case for any solution decaying faster than exponentially). The 3/2 exponent is independent 
of the space dimension; therefore the behaviour of the large c tail is singular for e — > 0, with an exponent jumping from 
3/2 for e > to 2 for e = (for dimensions larger than 1, the elastic system equilibrates and thus / is a Gaussian). 
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3. Numerical solution of the homogeneous Boltzmann equation 



The DSMC method gives access numericahy to the exact solution of the Bohzmann equation, and we restricted our 
study to the homogeneous situation. We have obtained the velocity distributions for various values of the restitution 
coefficient. Another powerful iterative method was recently introduced by Biben et al. Let us recall the idea of 
this method: 

The stationary velocity distribution can be obtained numerically through a direct iteration of equation (|^) . From an 
initial guess for the velocity distribution P{v) (a step function for example) the time evolution can be computed from 
equation (0) until the steady state is reached. Taking advantage of the v ^ —v symmetry of the velocity distribution, 
we only need to know the values of P for positive velocities. P{v) is then discretized from v = Q to v — [Ny — l)dv 
where typically = 1000 and dv = D^/tt / {n{l — a^))- The right hand side of equation can be estimated 

using Simpson integration, combined with a quadratic interpolation method to estimate the values of 



2z) - (1 - a)v' 



whose argument do not necessarily coincide with the velocity discretization. An implicit method is used to solve the 
diffusion equation: if p*+''* denotes the new estimate of P{i dv) at time t + dt, the left hand side of equation (|^) can 
be written in the time-velocity discretized space: 



pt+dt _ pt 



dt 



- D- 



t+dt 
i+l 



2K 



t+dt 



dv'^ 



pt+dt 



which leads to the following equation for P- 



t+dt. 



1 + 2D 



dt 

dv"^ 



pt+dt _ jj 



dt 
dv^ 



{Pitt + Pl-t) ^Pl+dt {RHS of (7) at site i and time t} . 



(13) 



We recognize a band tridiagonal matrix in the left hand side which can easily be inverted numerically to provide the 
new value of the velocity distribution at time t + dt . Normalization of the velocity distribution is enforced at each 
time step. 

Figure ^ shows a perfect agreement between the two methods, the iterative method allowing to reach a much higher 
precision (see the y-scale). The obtained distributions show strong deviations from the Gaussian, just as in the MD 
case. However, the distributions obtained by MD and DSMC agree only in the a ^ 1 limit, as was expected because 
of the spatial inhomogeneities appearing in the MD simulations. 




FIG. 2. Velocity distributions obtained by DSMC with 25000 particles (symbols) or by the iteration method (lines), in a 
log-linear scale, for restitution coefficients 0.1, 0.6, 0.9 and 0.99 (from top to bottom). Inset: same distributions on a linear 
scale, for restitutions 0.6 (stars) and 0.99 (squares). 

The measure of the fourth cumulant 02 (the first correction to the Gaussian) shows an excellent agreement between 
the DSMC data and the kinetic theory predictions from equation (^ij), on the whole range of inelasticities (see Fig. 

In the limit a ^ 1, 02 obtained in MD coincides with the prediction of Eq. ([TTI), namely, 02 —16/159. 
Moreover, the full velocity distribution function coincides with that obtained in DSMC (see Fig. H below). However, 
as a decreases, the MD results significantly deviate from their molecular chaos counterpart (see the inset of Fig. ||) 
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FIG. 3. Values of the second Sonine coefficient 02, obtained by measuring the fourth cumulant of /(c) in DSMC simulations 
(circles), together with the kinetic theory prediction Eq. |ll] (line). Inset: MD values for 02 (squares), and kinetic theory 
prediction (line). 

Moreover, plotting the velocity distribution versus either c^/^ or c"^, as in figure 0, shows that the high energy 
tail obtained by the DSMC method has a shape going from the predicted exp(— Ac /^) at large inelasticity to an 
exp{—Ac^) behaviour at low inelasticity. For intermediate values of a, a fit to a form exp(—Ac^) would lead to 
intermediate values of B. Since the real high energy limit has to follow an exp(— Ac'^/^) law M, this shows that for 
low £, this limit is far beyond reach of usual numerical methods, and emphasizes the fact that the range over which 
the large c limit is valid depends on the inelasticity. This is an important point since experiments have a limited 
precision, and the distribution function will have a practically vanishing weight much before this range is reached. In 
the next subsection we will see that the investigation of the a — > 1 limit and the use of the iterative method allows to 
understand the exp(— Ac"^) form obtained for a close to 1. 




c 



FIG. 4. Velocity distributions obtained by DSMC with 25000 particles, versus and in the inset versus c'^^^, for restitution 
coefficients 0, 0.1, 0.6, 0.707, 0.9, 0.99, 0.999 (from top to bottom): at low restitution coefficient /(c) shows an exp{~Ac''^^) 
behaviour, and an exp{—Ac^) behaviour as a goes to 1. 



B. Small inelasticity limit 



For small values of e = 1 — a, the Boltzmann equation takes the form 



dtP{v,t)^ DdlP{v,t) + ne / dw\v - w\P{w,t) 



P{v,t) + -{v-w)d^P{v,t) 



(14) 



The e — > limit can now be taken by introducing a scaled velocity x — {ne / D)^/^v. Using x, the Boltzmann 
equation yields an equation for the scaled distribution function 4>{x) = {ne/Dy/^Ps{v) (where Ps{v) is the stationary 
distribution limt^oo P^v, t)): 
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ay\x - y\(t>{y) 

— OO 

and we can integrate this equation twice to obtain 



(15) 



{x) = C exp 



Ay\y - (j){y) 



(16) 



Here C is a constant determined from the normalization condition J dx(l)[x) ~ 1. We have used the above equation 
to implement an iterative scheme to find numerically the corresponding velocity distribution. Moreover, as already 
pointed out in Q, one can easily see that the large |a;| limit is given by: 



(bix) = Ce-^l^l" , 
while at small x the function can be approximated by a Gaussian 



(t){x) = Ce-2 



(17) 



(18) 



with A = / da; |x|(/)(a;) « 0.785 determined from the numerical solution of Eq. ( |l6|) and C = C exp(— J dx\x\'^<j>{x)/6). 
The full numerical solution (displayed in Fig. ||) can also be investigated for locating the place where the crossover 



between the Gaussian and the exp (— |a:r/6) type behaviours takes place. We find that the crossover range deduced 
from comparing the asymptotics (RT() and dlq), Xa = 3A ~ 2.36, is actually in agreement with numerical observations 

of the full function. Thus returning to non-scaled velocities, we can see that the crossover velocity Wcr^ diverges in the 
e ^ limit as 



.,(1) 



D 



ne 



1/3 



(19) 



The important consequence of this result is that, for small dissipation, the cxp(— Aw"^) regime is reached for larger 
and larger velocities. 




FIG. 5. Symbols: velocity distributions obtained numerically at a = 0.999 by MD and DSMC simulations; the solid line is 
the numerical solution <f>{x) of Eq. (p^, corresponding to the a — > 1 limit, and here rescaled in order to have the same variance 
as the simulation data ((c^) = 1). The dotted line is the Gaussian distribution. 



At large but finite velocities, an effective exponent can be defined by 



d In V 



In < — In 



Piv) 



P(0) 



(20) 



corresponding to an apparent exp(— Aw") behaviour. The values of n for various inelasticities, obtained with the 
iterative method, are displayed in figure pi together with the a — > 1 limit. The effective exponent, starting from 2 at 



small velocities, increases and reaches a maximum at velocities that scale as e for small e (as Wcr )• The height of 
the maximum, ricr'' scales as 3 — ricr'' ^ e^^^. 
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One should note that the above scahng analysis explores only the v ^ e range of the velocity space. As already 
mentioned, for any fixed e, the large v limit of the distribution function is a stretched exponential 

P.(.)^exp|-^(^)'^'|.r|, (21) 

which is indeed observed numerically at relatively low values of a (see figure ^) . Comparing the arguments of the 
exponents in (|lj) and ( pT] ) one finds a crossover scale diverging as 

The effective exponent displayed in figure || indeed decreases at velocities larger than For large inelasticities, 
the n — > 3/2 limit of large velocities is observed; however, since vif' « Vcr^ /e^^'^, it becomes impossible to observe the 
asymptotics (^l]) for small inelasticities, even for almost realistic values of a like 0.95 (in experiments, a ~ 0.8 — 0.9). 

It has to be emphasized that this kind of behaviour is also observed in higher dimensions (for example, simulations 
of a two-dimensional heated granular gas with a = 0.8 yield an almost Gaussian velocity distribution, even if the 
predicted high energy behaviour is exp(— w^/^)); it is therefore to be kept in mind for the comparisons of models with 
experiments in which the available precision often does not allow to reach the predicted tails. 
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FIG. 6. Effective exponent defined by Eq. (|2(]|), for the solution of the Boltzmann equation obtained by the iterative method, 
and for the numerical solution of the a ^ I limit [upper dotted curve, obtained from the iterative resolution of Eq. (|l^)] . 

Finally, as already mentioned, it can be seen in Fig. |^ that for e — > 0, there is an excellent agreement between the 
single particle velocity distribution obtained in MD simulations (including both the space and velocity correlations), 
and that derived either from the asymptotic (e — > 0) distribution function (l){x) or from the Monte Carlo simulation 
of the Boltzmann equation. 

III. FREELY COOLING SYSTEM 
A. General considerations 

In the freely cooling system, no energy is injected and the temperature is monotonically decreasing with time. The 
first investigations of the one-dimensional freely evolving gas were undertaken in [^-^ ; it was shown by Molecular 
Dynamics simulations that, depending on the values of the number of particles and of the restitution coefficient, 
different instabilities could occur: e.g. at fixed number of particles N, if a is lower than a threshold, strong clustering 
occurs and leads to inelastic collapse [||. At larger a, the instability develops more slowly, and the inelastic collapse 
is avoided. The temperature then decays according to the rate equation dT/dt cx — T'^/^, i.e. T{t) '--^ however 
derived for an homogeneous system, whereas strong heterogeneities develop both in velocity and space coordinates; a 
wavy, oscillatory in time, perturbation appears in a "phase-space" plot (velocity versus position) with a tendency 
to form a bimodal velocity distribution. 
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The choice of a suitable quasi-elastic limit (where e ^ and A'^ cx 1/e to avoid the inelastic collapse) leads 
to a simplified Boltzmann equation [[tIJqHiiI, which can be understood using arguments from kinetic theory and 
hydrodynamics. This equation can be considered as formally exact as it has also been derived in from the exact 
BBGKY-like hierarchy governing the evolution of the distribution |^^. In this quasi-elastic limit, the velocity was 
observed to develop a two-hump structure reminiscent of the bimodal velocity distribution observed in Molecular 
Dynamics j?!^. Moreover, exact results were derived in the context of the above-mentioned limit, where it was shown 
in particular that to leading order in e, the velocity distribution consists of two symmetric Dirac peaks jl^ . 
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FIG. 7. Rescaled velocity distributions at large times at a = 0.85 obtained by MD with = 25000 (circles) and DSMC 
(squares) simulations. The solid line is the Gaussian distribution. In MD simulations, the inelastic collapse has been regularized 
by considering the same modified collision rule as in |l9| . 

Extensive MD simulations were carried out in using large sizes and probing large times, starting from an homo- 
geneous situation with an initial given velocity distribution. As long as the system is homogeneous, the temperature 
decays according to T{t) ~ t"^ As time evolves, space clustering of particles occurs and a decay is obtained 
Since the number of particles is large, inelastic collapse should then occur; this catastrophic event is avoided by 
imposing that the collisions between particles with relative velocity smaller than a given threshold are elastic (they 
can equally be chosen sticky), and it was checked that the results do not depend on the chosen threshold. In this 
respect, the authors of showed that the one-dimensional inelastic fluid belongs to the "universality class" of the 
sticky gas, and advocated a mapping to the Burgers equation to describe the appearing heterogeneities. Moreover, at 
large times the rescaled stationary velocity distribution /(c) was found very close to a Gaussian up to the available 
accuracy (see also Figure even if the mapping to the Burgers equation predicts an exp(— Ac'^) high energy tail. In 
fact, the bimodal structure of /(c) reported in |0Q can also be observed in this case during the transient homogeneous 
behaviour, during which spatial heterogeneities and correlations build up, as shown in |29| ; moreover, it can be clearly 
seen only by a convenient choice of the initial velocity distribution. The importance of spatial heteregeneities and 
correlations is emphasized in Fig. where the velocity distribution obtained following the prescription put forward 
by Ben-Naim et al. (that is essentially Gaussian) is compared to that obtained from the homogeneous Boltzmann 
equation pC| ]. The two-hump structure displayed by the latter appears for a > 0.8 and becomes more and more 
pronounced as a increases. In the next subsection, we will investigate in detail this structure in the e — > limit. 



B. Small inelasticity limit 



We have performed MD simulations using the two possibilities to avoid inelastic collapse mentioned in section III A 



Figure |[ left panel displays the results obtained for e — 5.10"'', at various stages of the evolution: at first the 
system remains homogeneous but the tendency to form a bimodal velocity distribution is rapidly obtained. The 
large time situation consists of two well defined clusters evolving with opposite velocities and yielding a sharply 
peaked bimodal velocity distribution (Figure ^, right panel). In this case, the overall kinetic energy E{t) ^ t^'^ 
decay consists in a series of plateaus since most of the dissipation occurs when the two clusters collide . 

On the other hand, with the regularization proposed in JTif , the duration of the transient behaviour increases 
with a, but the large time behaviour of the velocity distribution is always Gaussian. 
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FIG. 8. ie/i; Velocity-density scatter plots obtained in MD simulations with TV = 1000 and a = 0.9995. Each dot marks 
the location of a particle in the x — c plane, where x denotes the position in the simulation box and c denotes the velocity 
rescaled by the thermal velocity. Starting from an initial Gaussian distribution of velocities and random initial positions, the 
four snapshots correspond, from left to right and top to bottom, to four instantaneous configurations observed after respectively 
4.10^, 2.10*, 4.10^* and 6.10* collisions per particle. 

Right: Velocity statistics obtained in MD by averaging in the late time regime for the same initial situation and parameters as 
above, compared to its "mean-field" counterpart provided by DSMC simulations at the same inelasticity (a = 0.9995). 

It is striking to note that the two above procedures (small number of particles or elastic collisions at small enough 
velocities) lead to drastically different behaviours for the decay of the temperature, the spatial heterogeneities and the 
velocity distributions. It is moreover remarkable that the homogeneous solution of the Boltzmann equation captures 
the bimodality of /(c) (see Fig. H, right panel) associated with a strongly heterogeneous system. In order to gain 
insight into the approach to the limit e ^ 0, we therefore devote the remainder of this article to the analysis of the 
scaling properties of the homogeneous non-linear Boltzmann equation. We expect that for low enough inelasticity, 
/(c) tends towards two delta peaks at c = ±1, as predicted in |jl^. Performing DSMC simulations for smaller and 
smaller inelasticities allow to approach this regime, and Figure^ shows how the peaks become narrower as a increases. 




FIG. 9. Rescaled velocity distributions at q = 1 — e with e = 10 ,10 and 10 , obtained by DSMC simulations. 

As the system cools, the velocity distribution P{v,t) evolves into the Dirac distribution d{v). The above numerical 
results indicate that this distribution actually consists of two peaks located symmetrically around the velocity origin, 

for the rescaled velocity c 
id where the distributions 



at positions decaying like ±(et)~ . Moreover, it appears that the results displayed in Fig. p 
hide a self-similar structure, with a width of the peaks scaling like e^^"^, as evidenced in Fig, 
corresponding to different inelasticities collapse onto a master curve. The characteristic features of this master curve 
are investigated below and to this end, we return to the e expansion of the Boltzmann equation. (|lj), omitting the 
Fokker-Planck term Dd^P, since the fluid evolves freely in the present situation. 
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FIG. 10. Self-similarity of the rescaled distribution functions, for small inelasticities. 

From the evolution of temperature (T oc £^^<^^), it appears that a convenient scaling variable is a; = netv with a 
corresponding probability distribution function $ related to P{v,t) through P{v,t) — net^{nevt). To leading order 
in e, the Boltzmann equation is then brought into the simple form 



d(x$) 
da; 



= J dy\x-y\^{y) 



$(x) + -(x-y) — 



(23) 



the solution of which reads WM 



^{x)^ -[6{x-l)+5{x + l)]. 



(24) 



Looking for the self-similar structure of the peaks shown in Figures ^ and |l^ requires to push the e-expansion one 
order further compared to Eq. (|2S' 



= / dx'\x-x'\^x') 



d(a:$) 
dx 



and to consider solutions for $ of the form 



1 



L(l-£/2)2 

1 + b{e)x 



<& ( a; + ^ (a: - x')] - ^{x) 



1 - b{e)x 



(25) 



(26) 



In Eq. (|25|), the positive function -ip has the interpretation of the (e-rescaled) velocity distribution of left (or right) 
movers, and 6(e) = 6o + e''^! + £^'^^2- Substituting the scaling assumption for $ into ( |25| ) and performing the change 
of variables a; = — 1 + e'^y we obtain an non-linear intro-difFerential equation for ipiu)'- 



1 



a(e) [e'-%yi^)' ~ e^-^^')] - / dy'|y - v'\i^{y'my) + ^^iv - v'Wiv)) 



dy'\2~e%y + y'my') x 



(27) 



where terms of the form 'ip[—2e^'^ + y) have been neglected, anticipating that they will be exponentially small. Terms 
of order e^^^", e^+'' were equally neglected. Assuming again that t/i will have a sharp decay, we write that under the 
integrals |2 — e''{y + y')\ = 2 — e^iy + y'). Identifying on both sides of Eq. ( p7| ) terms of order e^"^"^ and e^~'^ leads 
respectively to 



(28) 



bo= I dj/V(y') 

which is the normalization condition, and 

i^'{y)(b,+ I dy'y'i^iy') 



(29) 
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which relates 61 to (y) (a constant function cannot be a solution). We choose to impose 61 = — (y)^ = 0, where the 
notation (...)^ stands for an average with weight function ^p. Then one notices that the expansion is consistent only 
to the condition that the e'^'^^ term can be cancelled by order e terms, which imposes that 



1 



(30) 



and we recover the exponent 1/3 that was needed to collapse the velocity distributions at several small inelasticities, 
as done empirically in Fig. ^ Finally we equate to the terms of order e to obtain the following integro-differential 
equation 



b24''{y) = / dyV(y') 



li^iy) i\y - y'\ -iy + y')) + \i^'{y) (|y - y'\{y - y') -{y + y'Y) + \r{y) 



which we integrate once, remembering that 60 = (1) 



ho^'{y) + 2b2^{y) + -V(2/) / dyV(y') (|y - y'|(y - y') - (y + y'f) = 0. 



(31) 



(32) 



We know from the direct analysis of the e — > limit of the scaling function that 60 = I7 which we use here on. At 
this stage we note that integrating the above equation and using that (y)^ = leads to 



262 = (y')v- 

Conversely, setting 62 = ^(y^)i/) will automatically enforce (y)^ = 0. We rewrite the equation for -0 in the form 

InV^(y) = InC - 262y + \ j dyV(y') [|y - yf - (y + y'f] . 

We now investigate the asymptotics of ip. The left tail of the distribution at large negative values of y reads: 

"1 



y 



ip{y)-C exp 



^y' + oii) 



(33) 



(34) 



(35) 



This sharp decay at y ^ —00 a posteriori justifies the approximations made in the course of the calculation. Note 
that omitting the first line in the rhs of Eq. ( p7| ) leads to exactly the same behaviour of the tail of the distribution. 
This has a physical interpretation: collisions between particles heading in the same direction can be neglected at large 
velocities. Opposite-velocity collisions are responsible for the form of the tail at large velocities. The o(l) represents 
contributions going exponentially fast to 0. 

The right tail y — > +00 decays exponentially fast as 



V'(y) ^ C exp {-2(y2)^ y + 0(1)} with C = C exp |-^(2/')vj 



(36) 



which again justifies the assumptions made so far. The o(l) again represents contributions going exponentially fast 
to 0. For completeness we mention the y ^ behaviour of the scaling function: 



0(y) = C'exp {-|(3s2 + \s\s)^ + y'{\s\)^) + 0{y')} , with C" = Cexp |i(|y|3 - y^) 



(37) 



As a side remark, the integro-differential equation for -0 can be cast in the form of an ordinary fourth-order 
differential equation 



(InV) 



(iv) 



(38) 



Finally, we note that numerical iteration of the integro-differential equation ( p^ ) converges extremely rapidly. This 
allows to determine the numerical constants C, (y^)^/i, (y^)^/i appearing in the asymptotics. The results obtained from 
this numerical scheme are compared with those of the DSMC method in Fig. and show a quantitative agreement 
which improves as e is lowered in DSMC, as expected (see the dotted curve at e = 10~^, closer to the asymptotic 
scaling form for ip than the dashed line corresponding to e = lO^'^). 
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FIG. 11. Comparison of the scaling function ip{y) (see text for definition) obtained within DSMC (e = 10 ^, dotted curve 



and e = 10 , dashed curve), with the solution of Eq. te2h corresponding to the quasi-elastic limit 



IV. CONCLUSION 

We have investigated the velocity statistics of one-dimensional granular fluids of inelastic particles, with a particular 
emphasis on scaling properties in the elastic limit, both in the absence of an external forcing and in a system heated 
by random accelerations. For the heated system, we showed that the expected high energy tail ^ exp(— Ac^/^) yields 
the correct asymptotic behaviour at finite inelasticity e, but this asymptotics is masked by a tail ^ exp(— Sc"^) for 
e ^ 0, with the rescaled crossover velocity between the two regimes scaling like e~^/^. This shows that the limits of 
high velocity and low inelasticity do not commute: if e — > is taken before the high energy limit, the distribution 
exhibits an asymptotic exp(— Ac'^) large c behaviour: 



, , c-»o< 

/(e,c) cx 
whereas lim fie.c) oc 



exp(-Ac^/^) for any e ^ 
exp(-74c^). 



(39) 
(40) 



Thanks to a high precision iterative scheme allowing to solve the homogeneous non-linear Boltzmann equation, we 
could obtain the velocity distribution over 30 orders of magnitude at arbitrary e, and thus define the apparent exponent 
n of the stretched exponential law for large c [/(c) cx exp(— Cc")]. However, even with such a precision, the crossover 
between the two behaviours (|3|) or (|o|) is difficult to investigate (see Fig. ||). 

For the freely evolving ID granular fluid, we have investigated in detail the structure and scaling properties of the 
two-hump velocity distribution exhibited by the solution of the homogeneous cooling state of the Boltzmann equation, 
both numerically and analytically. Such a bimodal distribution captures an essential feature of the velocity distri- 
butions obtained in Molecular Dynamics simulations for parameters hindering the inelastic collapse {i.e. extreme! 
small e or small systems). In this respect, a perturbative Sonine expansion in the spirit of that put forward in 
fails at small e, whereas such an expansion turned out to be remarkably accurate for the heated gas (see Fig. H). In 
both cases it would predict a non vanishing kurtosis correction a2 for £ — > 0, which is a peculiarity of d = 1; as soon 
as d > 1, a2 vanishes for small inelasticities, with or without forcing. 

It would be of interest to perform the same analysis for realistic space dimensions d > 1, and to quantify more 
precisely space and velocity correlations |5l|, for both the heated and unforced systems. 
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